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ABSTRACT 

We have investigated the role of the equation of state in resistive relativistic magne- 
tohydrodynamics using a newly developed resistive relativistic magnetohydrodynamic 
code. A number of numerical tests in one-dimension and multi-dimensions are carried 
out in order to check the robustness and accuracy of the new code. The code passes all 
the tests in situations involving both small and large uniform conductivities. Equations 
of state which closely approximate the single-component perfect relativistic gas are in- 
troduced. Results from selected numerical tests using different equations of state are 
compared. The main conclusion is that the choice of the equation of state as well as 
the value of the electric conductivity can result in considerable dynamical differences in 
simulations involving shocks, instabilities, and magnetic reconnection. 

Subject headings: magnetohydrodynamics (MHD) - methods: numerical - plasmas - 
relativistic processes 



1. Introduction 

Magnetic fields play an important role in determining the evolution of the matter in many 
astrophysical objects. In highly conducting plasma, the magnetic field can be amplified by gas 
contraction or shear motion. Even when the magnetic field is weak initially, the magnetic field can 
grow rapidly and influence the gas dynamics of the system. This is particularly important for the 
high-energy astrophysical phenomena related to strongly magnetized relativistic plasmas associated 
with objects such as active galactic nuclei (AGNs) (e.g., Urry &; Pavovani 1995), relativistic jets 
(e.g., Mirabel & Rodriguez 1999; Blandford 2002), pulsar winds (e.g., Gaensler & Slane 2006; 
Kirk et al. 2009), gamma-ray bursts (Zhang & Meszaros 2004; Piran 2005; Meszaros 2006), and 
magnet ars (e.g., Woods & Thompson 2006; Mereghetti 2008). 

The ideal magnetohydrodynamic (MHD) approximation is a good description of the global 
properties and dynamics of such systems well into their nonlinear regimes. In this limit the electri- 
cal resistivity rj = 1/a vanishes (infinite electrical conductivity). In this framework, many multi- 
dimensional ideal relativistic MHD (RMHD) codes have been developed to investigate relativistic 
astrophysical phenomena including fully non-linear regimes (e.g., Komissarov 1999; Koide et al. 
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1999; Komissarov 2001; Koldoba et al. 2002; Del Zanna et al. 2003; Leismann et al. 2005; Gammie 
et al. 2003; De Villiers & Hawley 2003; Anninos et al. 2005; Duez et al. 2005; Shibata & Sekiguchi 
2005; Anton et al. 2006; Mignone k Bodo 2006; Mizuno et al. 2006; Neilson et al. 2006; Del 
Zanna et al. 2007; Giacomazzo & Rezzolla 2007; Farris et al. 2008; Mignone et al. 2009; Beckwith 
& Stone 2011; Inoue et al. 2011). The ideal MHD limit provides a convenient form for solving 
the equations of RMHD and is also an excellent approximation for many relativistic astrophysical 
phenomena. However, in extreme cases such as binary mergers (the merger of two neutron stars 
or of a neutron star with a black hole) (e.g. Shibata <fe Taniguchi 2011; Faber & Rasio 2012) 
or the central engine of long GRBs (collapsar) (e.g., MacFadyen & Woosley 1999) the electrical 
conductivity can be small, and regions of high resistivity may appear. 

Quite often numerical simulations using ideal RMHD exhibit violent magnetic reconnection. 
The magnetic reconnection observed in ideal RMHD simulations is due to purely numerical resis- 
tivity, occurs as a result of truncation errors, and hence fully depends on details of the numerical 
scheme and resolution. Magnetic reconnection is one of the most important phenomena in astro- 
physics. It is highly dynamic, and it converts magnetic energy into fluid energy. The magnetic 
reconnection process has been invoked to explain flaring events (e.g., Lyutikov 2006; Giannios et al. 
2009) and magnetic annihilation (Coroniti 1990; Lyubarsky &: Kirk 2001) in relativistic plasmas. 
Therefore, numerical codes solving the resistive RMHD (RRMHD) equations and that allow control 
of magnetic reconnection according to a physical model of resistivity are highly desirable. 

Numerical simulation using the ideal RMHD equations is considerably easier than using the 
RRMHD equations because the equations become mixed hyperbolic with stiff relaxation terms. 
The pioneering work on resistive RMHD done by Komissarov (2007) solved the numerical flux by 
using the Harten-Lax-van Leer (HLL) approximate Riemann solver and by using Strang's splitting 
technique for the stiff relaxation terms. More recently, Palenzuela et al. (2009) have proposed 
a numerical method that solves the stiff relaxation terms in the equations by an implicit-explicit 
(IMEX) Runge-Kutta method, and Dionysopoulou et al. (2012) have extended the work of Palen- 
zuela et al. (2009) to 3D. A different approach has been taken by Dumbser Sz Zanotti (2009) 
who have applied the high order PnPm scheme to solving the resistive RMHD equations, and 
also Takamoto &: Inoue (2011) who have used the method of characteristics to solve the Maxwell 
equations accurately. Even more recently, a 3+1 resistive general relativistic MHD (GRMHD) code 
using mean-field dynamo closure has been developed by Bucciantini &: Del Zanna (2012). 

Plasma in the relativistic regime can have three major characteristics: the system has (1) 
relativistic fluid velocity (kinetic energy much greater than rest-mass energy), has (2) relativistic 
temperature (internal energy much greater than rest-mass energy), or has (3) relativistic Alfven 
speed (magnetic energy much greater than rest-mass energy). The second characteristic of relativis- 
tic temperature brings us to the issue of the equation of state (EoS) of the plasma. The EoS most 
commonly used in RMHD simulations is designed for plasmas with constant specific heat ratio (the 
so-called ideal EoS). However, this ideal EoS is valid only for plasmas with either ultra-relativistic 
temperature or non-relativistic temperature. The theory of relativistic perfect gases (Synge 1957) 
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has shown that the specific heat ratio cannot be constant if consistency with kinetic theory is re- 
quired. However, the exact EoS involves modified Bessel functions and is too complicated to be 
efficiently implemented in numerical codes. To get around this problem Mignone et al. (2005) in- 
troduced an approximate EoS given by a simple analytical formulation in the context of relativistic 
non-magnetized flows. This approximate EoS was applied in the context of relativistic MHD by 
Mignone & McKinney (2007). A different EoS approximation than that proposed by Mignone et 
al. (2005) has been proposed by Ryu et al. (2006). Clearly a determination of the effects of a 
difference in the choice of the approximate EoS is important to further advances in RRMHD. 

In this paper, we present the development of a new resistive RMHD simulation code including 
different realistic EoS approximations such as those proposed by Mignone et al. (2005) or by Ryu 
et al. (2006). This new RRMHD code is based on the ideal RMHD code RAISHIN (Mizuno et 
al. 2006; 2011) which uses a Godunov-type scheme to solve the conservation equations of ideal 
RMHD. In particular, we apply this new code to the role of the EoS in the resistive RMHD regime. 
We describe the basic equations of resistive RMHD in §2, three different equations of state are 
investigated in §3, and the numerical methods are described in §4. The various numerical tests in 
one-dimension and multi-dimensions are presented in §5. In §6 we conclude. 



2. Basic Equations of Resistive Relativistic MHD 

We have considered n a to be the time-like translational killing vector field in a flat (Minkowski) 
space-time, so n a = (—1,0,0,0), where we use Greek letters that take values from to 3 for the 
indices of 4D space-time tensors, while Roman letters take values from 1 to 3 for the indices of 3D 
spatial tensors. We use the speed of light c = 1 and Lorentz-Heaviside notation for electromagnetic 
quantities, so that all \/4vr factors disappear. 

The total energy-momentum tensor T a/3 describing a perfect fluid coupled to an electromag- 
netic field is defined as 

rpa/3 _ rpaj3 , rpa/3 /-. \ 

1 fluid ' EM' \ L ) 

The first term is due to matter: 

Tff uid = phu^+pg^, (2) 

where u a is the fluid four- velocity, while h (= 1 + e + p/p), p, p and e are the enthalpy, the proper 
rest mass density, the gas pressure, and the specific internal energy as measured in the fluid rest 
frame. The second term comes from the electromagnetic field: 

Tf M = F^FS-\(F^F, 1/ )g^, (3) 

where F a/3 , and its dual *F al3 are the Maxwell and Faraday tensors of the electromagnetic field 
given by 

*F al3 = n a B p - n !3 B a + n v e ua ^E^. (5) 



-4 - 



E a and B a are the electric and magnetic fields as measured by an observer moving along any- 
time-like vector n a , while e al3tJ,u = ^/—ge a p^ v is the Levi-Civita alternating tensor of space-time 
and € a pnv is the four-dimensional Levi-Civita symbol. 

In the global inertial frame with time-independent coordinate grid, the full system of Euler 
and Maxwell's equations are 



d t D + V • Dv = 0, 


(6) 


d t m + V • n = 0, 


(7) 


d t T + V • Y = 0, 


(8) 


3(£-VxB + V$ = -j, 


(9) 


d t B |Vx£ + V$ = 0, 


(10) 


d t ^ + V • E = q - k^, 


(11) 


<9 t $ + V • B = -k$, 


(12) 


d t q + V-j = 0, 


(13) 


q is the charge density, k is the damping 


rate parameter, and 



the conserved variables 

D = fry, (14) 

m = phj 2 v + E x B, (15) 

T = p h 7 2 -p+^(E 2 + B 2 ) (16) 

express the relativistic mass density, the momentum density, and the total energy density. Here, 
v is the velocity measured by an inertial observer and 7 = — v 2 is the Lorentz factor. The 
energy flux density and the momentum flux density can then be given by 



Y = phj 2 v + ExB, (17) 

9- (18) 



n = -EE - BB + ph-/ 2 vv + 



\(E 2 + B 2 ) +P 



An equation of state (EoS) is needed to close the system, we have adopted a variable EoS (e.g., 
Mignone et al. 2005; Mignone & McKinney 2007; Ryu et al. 2006). The details of the variable EoS 
are explained in next section. 

Eqs. (9)-(12) evolve the augmented Maxwell's equations which contain two additional fields ^ 
and <I> to control the system dynamics. In this approach, the two scalar fields ^ and $ indicate devi- 
ations of the divergence of the electric and magnetic fields from the values prescribed by Maxwell's 
equations, propagate at the speed of light, and decay exponentially over a time-scale ~ 1/n when 
the damping rate parameter k > 0. Following previous studies (Komissarov 2007; Palenzuela et 
al. 2009; Dumbser & Zanotti 2011), we have adopted the so-called hyperbolic divergence-cleaning 
approach used in the context of ideal MHD (Dedner et al. 2002). 
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The system of Eqs. (6)- (13) is closed by means of Ohm's law. Ohm's law for relativistic 
plasmas can be very complicated (e.g., Lichnerowicz 1967; Ardavan 1984; Blackman & Field 1993; 
Gedalin 1996; Melatos & Melrose 1996; Punsley 2001; Meier 2004). In this paper, we consider 
only the simplest kind of relativistic Ohm's law that assumes an isotropic plasma resistivity (e.g., 
Komissarov 2007; Palenzuela et al. 2009; Zenitani et al. 2010; Takamoto & Inoue 2011; Takahashi 
et al. 2011). In covariant form, the four-vector of the electric current is obtained from 

T a = a F a P U p + q u a , (19) 

where a = l/i] is the conductivity, r\ is the resistivity, and qo = —I a u a is the electric charge density 
as measured in the fluid flame (Lichnerowicz 1967; Blackmas k, Field 1993). In a special relativistic 
inertial frame, we find 

j = a-y[E + v x B - (E ■ v)v] + qv . (20) 
In the fluid rest frame, this equation becomes 

3 = ctE. (21) 

The ideal MHD limit of Ohm's law is given by the limit of infinite conductivity (a — > oo). In this 
limit Eq. (20) reduces to 

E + v x B- (E-v)v = 0. (22) 
Splitting this equation into components normal and parallel to the velocity vector give 

8 t E\\ + a-f[E\\ - (E ■ v)v] = 0, (23) 
d t E± + a-f[E± + v x B] = 0. (24) 

From these equations, one obtains the well-known ideal MHD condition 

E = -v x B. (25) 

In this limit the electric field is orthogonal to both magnetic and velocity fields. 



3. Equations of State 

An EoS relating thermodynamic quantities is required to close the system of eqs. (6)-(13). In 
general, an EoS is written as 

h = h(p,p), (26) 
and general forms for the polytropic index n and the sound speed c s are given by 

dh 9 p dh , „, 

n= i%- x '*—kw (27) 
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The most commonly used EoS, a constant T-law (ideal) EoS, is given by 

h = l + y^—Q, (28) 

where T is the constant specific heat ratio and 9 = p/p is the temperature. The sound speed is 
calculated from 

4^. (29) 

The constant T-law EoS may be applied correctly to a plasma with non-relativistic temperature 
where T = 5/3 or to a plasma with an ultra-relativistic temperature where V = 4/3. However, in 
the high-temperature limit, i.e., — > oo with T > 4/3, the sound speed exceeds relativistic limit 
(c s > l/\/3). Moreover a constant T-law EoS is not consistent with relativistic kinetic theory, the 
so-called Taub's fundamental inequality, which requires the specific enthalpy to satisfy 

(h-Q)(h-4&) > 1. (30) 

This rules out a constant T-law EoS with T > 4/3, if applied to < 6 < oo. 

The theory of single-component perfect gases in the relativistic regime shows that the specific 
enthalpy is a function of the temperature © = p/p only, and has the form (Synge 1957) 

h - K3il/Q) (31) 

where K 2 and K% are the 2 nd and 3 rd order modified Bessel functions of the second kind respectively. 
Using an equivalent T eq = [h — l)/(h — 1 — 0) in the non-relativistic temperature limit (G — > 0) 
yields T eq — > 5/3, and in the ultra-relativistic temperature limit (O — > oo) yields T eq — > 4/3 (see 
Fig. la). However, this EoS requires extra computational costs because the thermodynamics of 
the fluid is expressed in terms of the modified Bessel functions (Falle & Komissarov 1996). 

Recently, Mignone et al. (2005) proposed an EoS, the so-called TM EoS, that follows eq. (31) 
well. The TM EoS, which was first introduced by Mathews (1971), is given by 



and the sound speed is calculated from 



2 6 5/i - 80 



The TM EoS corresponds to the lower bound of Taub's fundamental inequality, i.e., (h — Q)(h — 
46) = 1, and produces the correct asymptotic values for T eq . 

Ryu et al. (2006) proposed an EoS which is a simpler algebraic function of 6, hereafter referred 
to as the RC EoS, that satisfies Taub's inequality for all 0. The RC EoS is given by 

p 3p + 2p 602 + 40 + 1 

or h = 2 — — , (34) 



e-p 9p + 3p 30 + 2 
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and the sound speed is calculated from 




see 



Fig. 1. — Equivalent T (left), specific enthalpy (middle), and sound speed (right) as functions of 
the temperature Q = p/p. Different lines correspond to: constant T-law EoS with T = 5/3 (dotted 
lines), constant T-law EoS with T = 4/3 (dashed lines), TM EoS (dash-dotted lines) and RC EoS 
(dash-two dotted lines). For comparison Synge's EoS has been plotted as the solid lines. 

Figure 1 shows the equivalent T, the specific enthalpy and the sound speed as a function of 
for TM EoS, RC EoS, Synge's EoS as well as a constant T-law EoS with T = 5/3 and 4/3. 
The specific enthalpy and the sound speed computed using the TM EoS and the RC EoS are well 
matched to Synge's EoS and cannot be distinguished on the plots. The approximations to Synge's 
EoS such as the TM EoS and RC EoS are hereafter referred to as approximate EoSs. 



4. Numerical Method 

A well-known and challenging feature of the system of equations (6)- (13) is that they have 
source terms for the evolution of the electric field that become stiff in the high conductivity (low 
resistivity) limit. Following Komissarov (2007), we will use the Strang-splitting technique (Strang 
1968). 

The system of equations (6)-(13) can be written as a single phase vector equation 

dU(P) dF m {P) 
dt + dx m ~ 



(36) 



where 
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are the vectors of conserved variables, primitive variables and sources, respectively, and 



fTfV 
ym 

-e imk B k + ^g im 
3 imk E k + $g im 

E m 



\ 



V 



is the vector of numerical fluxes, where e y ' fe is the Levi-Civita alternating tensor of space. 



The source term can be split into two parts 



Sa(P) 





0* 


—qv l 
0* 
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and S b (P) 



J 



V 






-fc 

o i 





where 



j c = aj[E + v x B - (E- v)v\ (37) 

is the conductivity current. The source term Sb is a stiff relaxation term that requires special 
care to capture the dynamics in a stable and accurate manner. In the Strang time-step splitting 
technique, firstly the solution is advanced using the stiff-part equations 

dU{P) 



dt 



S b (P) 



(38) 



over the half time-step, At/2. Secondly, advance of the non-stiff part of the equations is made via 
second-order accurate numerical integration of 



au ( p) + ™ = 



dt 



Q x m 



(39) 
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over the full time step. Thirdly, again the solution is advanced by the stiff-part equations over the 
half-time step. 

Time advance of the non-stiff part of the equations is given by 

TT tt | A + -^-1/2^+1/2 ~ i ? m+l/2,n+l/2 ( . , 

U n+1 - U n + At 2^ 1- Atb a:n+1 /2, (40) 

m=l 

where U n represents the cell-centered conserved variables at t = t n , U n+ \ represents the cell-centered 
conserved variables at t = t n + At, S a ^ n+ i/2 represents the cell-centered non-stiff source term at 
t = t n + At/2, F m+1 / 2t n + i/2 is the numerical flux though the right-hand side cell interface and 
^m-i/2,n+i/2 is the numerical flux though the left-hand side cell interface in the direction of x m at 
time t = t n + At/2, Ax m is the cell size in this direction, and rid is the number of spatial dimensions. 
The non-stiff sources and numerical fluxes at the half-time step are calculated from 

-v.+%£ Fm - i/2 -;;f" +i/ ' 2 ''' + < 4i > 

m=l 

The numerical fluxes at the cell-interface F m+1 / 2n are calculated using the simplified Harten-Lax- 
van Leer (HLL) approximate Riemann solver (Harten et al. 1983, Komissarov 2007) where the 
maximum characteristic speed of the system in each direction equals the speed of light. The left- 
hand and right-hand states of each cell interface using the HLL approximate Riemann solver are 
computed from various reconstruction schemes as in our ideal RMHD code (Mizuno et al. 2006; 
2011). In this paper, we use a piecewise linear method (PLM) reconstruction such as the minmod 
slope-limited linear interpolation scheme or the Monotonized Central (MC) slope-limited linear 
interpolation scheme as these are the simplest reconstruction schemes that capture a shock sharply. 
The resulting scheme is second-order accurate in time and space. 

Following the work of Komissarov (2007), the split evolution equations of the electric field, 
eqs. (23) & (24), can be solved analytically 



E n = £f?exp 



7 . 



(42) 



E ± = E* ± + (Efi ± -E* ± )exp[-ajt], (43) 

where E* ± = —v x B and suffix indicates the initial component. The stiff-part equations related 
to the two scalar fields ^ and <I> are also solved analytically with solution 

^ = * exp[-Kt], (44) 
$ = <& exp[-/ci], (45) 

where and are the initial values of ^ and <!>. 

In order to evolve this system of equations, the numerical fluxes F m must be computed at 
each time-step. These fluxes depend on the primitive variables P, which must be recovered from 
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the evolved conserved variables U. In conserved variables, E and B can be calculated by evolving 
Maxwell's equations. However, it is more stable to evolve the stiff part equations (42) - (43) during 
the primitive recovery process when a becomes large (Palenzuela et al. 2009). The primitive 
recovery procedure adapted to an approximate EoS such as the TM EoS and the RC EoS follows 
that used by Palenzuela et al. (2009). 



5. Numerical Tests 

In this section, one-dimensional and two-dimensional tests are presented. Three one-dimensional 
tests have been used to validate our new resistive relativistic MHD code in different regimes. A 
one-dimensional shock-tube test and three two-dimensional tests have been used to investigate the 
effect of different EoSs. The damping coefficient of the hyperbolic divergence cleaning is set to 
k = 1. The magnetic field is divergence-free and charge is preserved at the truncation error level. 



5.1. One-dimensional tests 

5.1.1. Large amplitude CP Alfven wave test 

This test consists of the propagation of a large amplitude circularly-polarized Alfven wave 
along a uniform back-ground magnetic field Bq in a domain with periodic boundary conditions. 
The exact solution is given by Del Zanna et al. (2007) in the ideal MHD limit, and was used as an 
ideal-MHD limit test problem by Palenzuela et al. (2009) and Takamoto & Inoue (2011). Here we 
use conditions similar to previous studies with 

{B y , B z ) = ( A B {cos[k(x - v A t],sm[k(x - v A )t}), (46) 

(v y ,v z ) = —j^-(B y , B z ), (47) 
Bo 

where B x = Bq, v x = 0, k is the wave number and ( A is the amplitude of the wave. The special 
relativistic Alfven speed v A is given by 

^-mlMW'-Gn^^r (48) 

For this test we use initial parameters p = p = 1 and Bq = 0.46188, the Alfven velocity 
v A = 0.25c, and we adopt a constant gamma-law EoS with T = 2. Following Palenzuela et al. 
(2009), we use a high uniform conductivity a = 10 5 with three different resolutions of 50, 100 and 
200 cells covering the computational domain x G [—0.5,0.5]. Figure 2 shows the numerical results 
at t = 4 (one Alfven crossing time) for the three different resolutions. This result shows that the 
new resistive RMHD code reproduces ideal relativistic MHD solutions when the conductivity a is 
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Fig. 2. — Magnetic field component B y in a large- amplitude CP Alfven wave test using three 
different resolutions iV = 50 (dotted), 100 (dashed) and 200 (dash-dotted) at t = 4. The solid line 
shows the exact solution. The numerical results are in good agreement with the analytical one (the 
highest resolution is excellent). 

high. The L\ norm errors of the magnetic field component B y in this test are shown in Figure 3. 
The numerical result is slightly shallower than 2 nd order convergence. This is likely caused by the 
periodic boundary. 

0.100 



o.oro 



0.001 

w too 1000 

Fig. 3. — L\ norm errors of the magnetic field component B y in a large-amplitude CP Alfven wave 
test using the three different resolutions N = 50, 100 and 200. 

5.1.2. ID self- similar current sheet test 

This test problem has been used for moderate resistivity cases (e.g., Komissarov 2007; Palen- 
zuela et al. 2009; Takamoto & Inoue 2011). In this test problem the magnetic pressure is 




2nd order 
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much smaller than the gas pressure everywhere. The magnetic field configuration is given by 
B = [0, B y {x, t),0], where B y (x,t) changes sign within a thin current layer of thickness Al. An 
initial solution is provided in equilibrium with p = constant. The evolution of this thin current 
layer is a slow diffusive expansion due to the resistivity and described by the diffusion equation 



d,B„ 



-d 2 x B y 



0. 



(49) 



As the thickness of the layer becomes much larger than Al the expansion becomes self-similar with 



B !l {.r.l)=B () ~J~ 



(50) 



where x = t/x 2 and erf is the error function. This analytic result can be used for testing the 
moderate resistivity regime. 

In the test problem, we have chosen an initial solution at t = 1 with p = 50, p = l,E = v = 
and a = 100 (ry = 1/er = 0.01). A constant gamma-law EoS with T = 2 is used. The computational 
domain is uniform with 200 cells in [—1.5, 1.5]. The numerical simulation is evolved up to t = 10 and 



0.5 - 




Fig. 4. — Magnetic field component B y in a Self-similar current sheet test. The dotted and dashed 
lines are indicated the analytical solution at t=l and t=10. The solid line shows the numerical 
solution at t = 10. The numerical solution is in excellent agreement with the analytical one. 

then the numerical solution is compared in Figure 4 with the analytical solution. The numerical 
and analytical solutions cannot be distinguished on the plot. This indicates that the moderate 
resistivity regime is well described by the code. 



5.1.3. ID shock-tube tests 



As a first shock-tube test in the restive relativistic MHD regime, we consider a simple MHD 
version of the Brio and Wu test as in Palenzuela et al. (2009) and Takamoto & Inoue (2011). The 
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initial left and right states are separated at x = 0.5 and are given by 

(p L ,p L ,B^) = (1.0,1.0,0.5) (51) 
(p R ,p R ,B R ) = (0.125,0.1,-0.5). (52) 

All other fields are set to 0. We use a constant T-law EoS with T = 2. The computational domain 
covers the region x £ [—0.5, 0.5] with 200 cells. 

Figure 5 shows the numerical results at t = 0.4 for conductivities a = 0, 10, 10 2 , 10 3 , 10 5 . 
The exact solution to the ideal RMHD Riemann problem was found by Giacomazzo & Rezzolla 
(2006). When B x = 0, the solution contains only two fast waves, a left-moving rarefaction wave 
and a right-moving shock with a tangential discontinuity between them. The results show that the 




-0.4 -0.2 0.0 0.2 0.4 




-1.01 , , , , . , , , . , . , . , , . . ,_ 
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Fig. 5. — (a) Density and (b) magnetic field component B y in the simplified Brio & Wu shock-tube 
test. Different lines indicate different conductivities: a = (orange solid), 10 (green dash-two- 
dotted), 10 2 (red dash-dotted), 10 3 (purple dashed), 10 5 (blue dotted). The black solid line shows 
the exact solution in the ideal RMHD case. 

solution smoothly changes from a wave-like solution for a = to the ideal-MHD solution for high 
conductivity a = 10 5 . Note that for a = the solution describes a discontinuity propagating at 
the speed of light corresponding to Maxwell's equations in vacuum. This result is nearly the same 
as test results for other codes (Palenzuela et al. 2009; Takamoto &; Inoue 2011). 
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Palenzuela et al. (2009) reported that Strang's splitting technique became unstable for mod- 
erately high conductivity in this shock tube test, and they suggested using an implicit method. 
However, Takamoto & Inoue (2011) found that this instability is not related to Strang's splitting 
technique but instead to the calculation of the electric field during the primitive recovery proce- 
dure. In this new resistive RMHD code, the shock tube test is solved stably using Strang's splitting 
technique even when a ~ 10 6 . 

Balsara Test 2 (Balsara 2001) is used as a second shock tube test to investigate the effect of 
different EoSs in the resistive relativistic MHD regime. In ideal RMHD, Mignone & McKinney 
(2007) have already performed this test to check the effect of different EoSs. In this test the initial 
left and right states are separated at x = 0.5 and are given by 

(p L ,p L ,B^B^B^) = (1.0,30.0,5.0,6.0,6.0) (53) 
(p R ,p R ,B R ,B R ,B R ) = (1.0,1.0,5.0,0.7,0.7). (54) 

The computational domain covers the region x S [—0.5, 0.5] with 800 cells. This test shows that a 
mildly relativistic blast wave propagates to the right with maximum Lorentz factor of 1.3 < 7 < 1.4. 

The numerical results at t = 0.4 for the constant T-law EoS with T = 5/3, the TM EoS and the 
RC EoS using conductivities a = 0, 10, 10 2 , 10 3 are shown in Figures 6, 7 and 8, respectively. The 



Gamma-law EoS 




Fig. 6. — (a) Density, (b) gas pressure, (c) velocity component v x , (d) velocity component v y , (e) 
magnetic field component B y and (f) Lorentz factor in the Blasara Test 2 (mildly relativistic blast 
wave) at t = 0.4 using an ideal EoS with T = 5/3. Different lines indicate different conductivity: 
a = (purple dash- two-dotted), 10 (green dash-dotted), 10 2 (red dashed), 10 3 (blue dotted). The 
black solid line shows the exact solution in the ideal RMHD case. 

solutions show fast and slow rarefaction waves, a contact discontinuity, and slow and fast shocks 
from left to right. In these cases, no rotational discontinuity is seen. The results obtained from the 
TM EoS and RC EoS cases are considerably different from the constant T-law EoS with T = 5/3. 
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Fig. 8. — The same as in Fig. 6 but using the RC EoS. 



In the approximate EoS cases, the rarefaction waves and shocks propagate with smaller velocities. 
This is predicted from the lower sound speed in the approximate EoS cases relative to overestimated 
sound speed in the ideal EoS case with T = 5/3. Behind the slow shock, the approximate EoS cases 
have a higher peak density, which follows from the previous considerations. These properties are 
consistent with those in the ideal RMHD case (Mignone & McKinney 2007). On the other hand, 
the results obtained from TM EoS and RC EoS cases are very similar at our numerical resolution. 
This similarity reflects the similarity in the distributions of specific enthalpy (see Fig. 1). Again the 
results show a smooth change from a wave-like solution for a = towards an ideal-MHD solution 
for the highest conductivity, a = 10 3 , for all the different EoS cases. The differences between the 
approximate EoS cases and the constant T-law case are larger at higher a where the approximate 
EoS cases approach the ideal MHD case more slowly. Even for low conductivity, i.e, a = 10, we 
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clearly see a difference between the ideal EoS case with T = 5/3 and the approximate EoS cases. 

5.2. Two-dimensional tests 

5.2.1. The Cylindrical Explosion 

We now consider tests involving shocks in multi-dimensions. Firstly we choose a test involving 
a cylindrical blast wave expanding into an initially uniform magnetic field. This is a standard 
test for ideal relativistic MHD codes even though there is no exact solution because this test will 
reveal subtle bugs and potential weaknesses in the numerical implementation. For this test, we use 
a Cartesian computational domain (x, y) £ [—6, 6] with 200 uniform cells in each direction. The 
initial explosion is initialized by setting the gas pressure and density to p = 1 and p = 0.01 within 
a cylinder of radius r < 0.8 centered on the origin. In an intermediate region 0.8 < r < 1.0, the 
pressure and density decrease exponentially to that of the ambient gas which has p = p = 0.001. 
The initial magnetic field is uniform in the x-direction with B = (0.05,0,0). The other quantities 
are set to zero (i.e., v = E = q = 0). 

Figure 9 shows the magnetic field components B x and B y at t = 4 using the ideal EoS with 
T = 4/3 and the TM EoS. The ideal-MHD simulation is performed using a high conductivity of 
a = 10 5 . The results are qualitatively similar to those reported in previous studies in ideal RMHD 
(Komissarov 1999; Leismann et al. 2005; Neilsen et al. 2006; Del Zanna et al. 2007; Mizuno et 
al. 2011) and RRMHD (Komissarov 2007; Palenzuela et al. 2009). The results obtained from the 
ideal EoS with T = 4/3 and the TM EoS are qualitatively very similar. This means that an ideal 
EoS with T = 4/3 satisfactorily captures all the shock properties. 

Figure 10 shows one-dimensional profiles of the gas and magnetic pressure along the y-axis for 
conductivities of a = 0, 10, 10 2 , 10 3 , 10 5 using the ideal EoS with T = 4/3 and the TM EoS. At high 
conductivity a > 10 3 , we do not see any significant difference. This means that high conductivity 
recovers the ideal-MHD solution. As the conductivity decreases, the maximum gas and magnetic 
pressure decrease. Of course there is no magnetic pressure increase for a = 0. 

5.2.2. Kelvin- Helmholtz instability test 

We present calculations of the linear and nonlinear growth of the two-dimensional Kelvin- 
Helmhotz instability (KHI) to investigate the effect of conductivity and the EoS on the development 
of turbulence in the resistive relativistic MHD regime. 

Initial conditions for this test are taken from a combination of previous studies in ideal RMHD 
(Bucciantini & Del Zanna 2006; Mignone et al. 2009; Beckwith & Stone 2011). The shear velocity 
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(a) Gamma-law EoS (b) 




(c) TMEoS (d) 




Fig. 9. — Magnetic field components B x (left panels) and B y (right panels) for the cylindrical 
explosion test at t = 4 using a ideal EoS with T = 4/3 (upper panels) and TM EoS (lower panels). 



profile is given by 

v th tenh ify>0.0 

V X = < V , „ r \ (55) 

\ -^tanh(^) ify<0.0 

Here, a = 0.01 is the characteristic thickness of the shear layer, and v s h = 0.5 corresponds to a 
relative Lorentz factor of 2.29. The initial uniform pressure is p = 1.0. The density is initialized 
using the shear velocity profile, with p = 1.0 in regions with v s h = 0.5 and p = 10 -2 in regions 
with v s h = —0.5. The magnetic field components are given in terms of the poloidal and toroidal 
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Fig. 10. — Gas pressure P gas (left panels) and magnetic pressure P mag (right panels) for the 
cylindrical explosion test at t = 4 using an ideal EoS with T = 4/3 (upper panels) and a TM EoS 
(lower panels). The different lines show conductivity cases: a = (purple dash- two-dotted), 10 
(green dash-dotted), 10 2 (red dashed), 10 3 (blue dotted), and 10 5 (black solid). 

magnetization parameters fi p and fit as 



{B x , B y , B z ) = (^2AW 0, y%Hp), 



(56) 
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with fi p = 0.01 and fit = 1.0. The instability is seeded by a single mode perturbation of the form 



Aov s h sin(27rx) exp 
— Aov s h sin(27rx) exp 



y-o.5 



y+0.5 



if y > 0.0 
if y < 0.0 



(57) 



Here, Aq = 0.1 is the perturbation amplitude and a = 0.1 is the characteristic length scale over 
which the perturbation amplitude decreases exponentially. The computational domain covers x G 
[-0.5,0.5], y G [-1,1] with 256 x 512 cells. 
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Figure 11 shows the perturbation amplitude (Av y 
poloidal magnetic field (B po i = + £? 2 ) as a function of time for conductivities of a 

10 2 , 10 3 , 10 5 using an ideal EoS with T = 4/3 and using the TM EoS. All cases show an initial 
linear growth phase. Except for a = 0, the ideal EoS and the TM EoS have almost the same 
growth rate and the maximum amplitude is reached at t ~ 2. The maximum amplitude indicates 
the transition from the linear to the nonlinear phase. The a = cases exhibit a lower growth rate 
and later transition to the nonlinear phase than the higher conductivity cases. Thus, differences in 
the conductivity and EoS do not affect the growth of KHI, except for a = 0. 
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Fig. 11. — Evolution of the amplitude of the perturbation (upper panels) and volume-averaged 
poloidal field {B po i) (lower panels) as a function of time for the Kelvin- Helmholtz instability test 
using a ideal EoS with V = 4/3 (left panels) and TM EoS (right panels). The different lines 
indicate different conductivity cases: a = (purple dash- two-dotted) , 10 (green dash-dotted), 10 2 
(red dashed), 10 3 (blue dotted), 10 5 (black solid). 



Poloidal field amplification via stretching due to the main vortex developed by KHI follows the 
growth of KHI (see Fig. 12). In high conductivity cases, poloidal field amplification is very large, 
an increase by almost one-order of magnitude. Saturation in the poloidal field amplitude occurs 



-20- 



latter than the transition from the linear to the non-linear KHI growth phase. This means that 
magnetic field amplification via stretching continues even after KHI is fully developed. Poloidal 
field amplification is weaker and saturation occurs earlier when the conductivity is low. Larger 
poloidal field amplification occurs for the TM EoS case than for the ideal EoS case. Therefore we 
find that magnetic field amplification via stretching due to the main vortex developed by KHI is 
strongly affected by the conductivity and the EoS. 

Figures 12 and 13 sh ow the time evolution of the density and the poloidal to toroidal field ratio 
(B pol /B tor = yjBl + B$/B z ) for high conductivity a = 10 5 , using the ideal EoS with T = 4/3 
(Fig. 12) and using the TM EoS (Fig. 13). Both cases show formation of a main vortex by growth 
of KHI in the linear growth phase. In the ideal EoS case, a secondary vortex appears, although 
not fully developed. However, development of a secondary vortex is not found in the TM EoS case. 
Beckwith & Stone (2011) found that a secondary vortex did not appear even in very high resolution 
simulations using the HLL approximate Reimann solver in ideal RMHD. In our simulations, we also 
used the HLL approximate Riemann solver to calculate the numerical flux, but do find a secondary 
vortex in the ideal EoS case. The difference is likely the result of the reconstruction scheme used 
here and the different reconstruction scheme used by Beckwith &; Stone. In the nonlinear phase the 
main vortex is distorted and stretched. The magnetic field is strongly amplified by shear motion in 
the vortex in the linear phase and by stretching in the nonlinear phase. As the mixing layer grows 
the field lines are bunched into a filamentary like stretched structure. In the TM EoS case, the 
vortex becomes strongly elongated along the flow direction. The structure created in the nonlinear 
phase is very different in the ideal EoS and the TM EoS cases. 

The field amplification structure for different conductivities from a = to 10 3 is shown in 
Figure 14. As seen in the time evolution of the averaged poloidal field shown in Fig. 11, the 
magnetic field amplification is weaker when the conductivity is low. Field amplification is a result 
of fluid motion in the vortex. In the high conductivity case, the magnetic field follows the fluid 
motion, like ideal MHD, and is strongly twisted. When the conductivity declines, the magnetic 
field is no longer strongly coupled to the fluid motion. Therefore the magnetic field is not strongly 
twisted. In the case using the TM EoS, we see the same trend for different conductivity and do 
not show the result here. 



5.2.3. Relativistic Magnetic Reconnection test 

The final test involves relativistic magnetic reconnection. Pioneering work on relativistic mag- 
netic reconnection using a resistive relativistic MHD code and 2D simulations was performed by 
Watanabe & Yokoyama (2006) who considered Petschek-type reconnection. Zenitani et al. (2010) 
also have studied details of Petschek-type reconnection via 2D resistive relativistic MHD simula- 
tions. Zanotti & Dumbser (2011) have investigated the dependence of Petschek-type relativistic 
magnetic reconnection by performing 2D and 3D simulations over a broad range of conductivities 
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Fig. 12. — Density (upper panels) and the poloidal to toroidal field ratio (B po [/ B tor ; lower panels) 
for the Kelvin-Helmholtz instability test at t = 3, 7, & 10 for a = 10 5 using the ideal EoS with 
E = 4/3. White lines indicate magnetic field lines and the arrows show velocity vectors. 



and magnetizations. Takahashi et al. (2011) have studied Sweet-Parker type relativistic magnetic 
reconnection using 2D resistive relativistic MHD simulations. In this test, we present simulations 
of Petschek-type reconnection and investigate the effect of the EoS. 

We use initial conditions similar to that used in previous work (Watanabe & Yokoyama 2006; 
Zenitani et al. 2010; Zanotti & Dumbser 2011). The density and gas pressure are given by 



p = pb + Pm cosh 2 (y) 



(58) 
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Fig. 13.— The same as Fig. 12 but using the TM EoS. 



P = P6 + /x m cosh 2 (y), 



(59) 



where pb = Pb = 0.1 are the uniform density and gas pressure outside the current sheet, and 
[i m = -Bq/(27q) = 1.0 is the magnetization parameter. The velocity field is initially zero, hence 
7o = 1. The magnetic field changes orientation across the current sheet according to 

B x = B Q twh(y), (60) 

where -Bo is calculated from the magnetization parameter. The current distribution is given by 

j z = 5 cosh- 2 (y). (61) 
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Fig. 14. — The poloidal to toroidal field ratio (B po i/Bt or ; lower panels) for the Kelvin- Helmholtz 
instability test at t = 3 for (a) a = 0, (b) a = 10, (c) a = 10 2 , and (d) a = 10 3 using the ideal EoS 
with T = 4/3. 



Over the whole computational domain there is a small background uniform resistivity % = 1 /erj, = 
10~ 3 , except within a circle of radius r„ = 0.8 which defines a region of anomalous resistivity with 
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(62) 



where r = y 1 x 2 + y 2 . The electric field is calculated from the resistivity distribution as 



E z = VJz- 



(63) 



The computational domain is x G [—50,50], y G [—20,20] with 2000 x 800 cells, and outflow 
boundary conditions are used in both directions. 

Figure 15 shows the density, the x-component of the 4- velocity, "fv x , and the out-of-plane 
current j z at t = 100 using an ideal EoS with T = 4/3 and the TM EoS. In this figure we 
confirm the essential features of Petschek-type relativistic reconnection reported in previous work 
(Watanabe & Yokoyama 2006; Zenitani et al. 2010; Zanotti & Dumbser 2011). In both the ideal 
and TM EoS cases, we see similar time evolution and morphology. After an initial adjustment stage 
of t < 10, the reconnection process starts around the point at (x, y) = (0, 0) triggered by anomalous 
resistivity. The magnetic field shows a typical X-type topology. As a result of reconnection, the 
magnetic energy is converted into both thermal and kinetic energy. Two magnetic islands (so-called 
plasmoids) which correspond to the high-density region in Fig. 15 move in opposite directions (the 
figure shows only half of the simulation region and only one magnetic island) and are accelerated 
along the direction of the magnetic field. A fast reconnection jet is formed inside a narrow nozzle 
within a pair of slow shocks (Petschek slow shock). The reconnection jet collides with a plasmoid 
in front of the current sheet further downstream. The plasmoid is surrounded by strong currents 
(see Fig. 15e,f), which also correspond to the slow shocks (Ugai 1995). These slow shocks surround 
the plasmoid and are connected to the Petschek slow shocks. In the TM EoS cases, the plasmoid 
has a faster speed and than in the ideal EoS case with V = 4/3, and the plasmoid in the TM EoS 
case propagates further. 

The time evolution of the maximum outflow velocity (v Xtmax ), the volume-averaged magnetic 
energy (B 2 ) and the normalized reconnection rate is shown in Figure 16. The outflow gradually 
accelerates and nearly saturates by t ~ 60 with v x ~ 0.8c. The figure clearly shows that the 
outflow speed in the TM EoS case is slightly faster than in the ideal EoS case. The magnetic 
energy gradually decreases with time. Dissipated magnetic energy results in an increase to both 
the thermal and kinetic energy. Increase in the kinetic energy accompanies acceleration of the 
outflow (plasmoid). The normalized reconnection rate is defined as 1Z = E* /vA,i n Bi n ~ Vi n /v ou t, 
where E* is the electric field at the reconnection point and the upstream properties with subscript 
in are evaluated at (x, y) = (0, 3) (Zenitani et al. 2010)0. The reconnection rate saturates at about 
t = 50 in both cases. However, the TM EoS case has a larger reconnection rate than the ideal EoS 
case (1Z ~ 0.16 in the ideal EoS case with T = 4/3 and 1Z ~ 0.17 in the TM EoS case). Therefore 
the different EoSs lead to a quantitative difference in relativistic magnetic reconnection. 



Zanotti & Dumser (2011) and Takahashi et al. (2011) have used a different definition for late reconnection. 
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Fig. 15. — Density (upper panels), the ^-component of the four- velocity ^v x (middle panels), and 
the out-of-plane current j z (lower panels) in the relativistic magnetic reconnection test at t = 100 
using the ideal EoS with T = 4/3 (left panels) and the TM EoS (right panels). White lines indicate 
magnetic field lines and the arrows show velocity vectors. 



6. Summary and Conclusions 

The role of the EoS in resistive relativistic MHD using a newly developed resistive relativistic 
MHD code has been investigated. A number of numerical tests in ID and multi-dimensions have 
been performed to check the robustness and accuracy of the new code. All of the tests show the 
effectiveness of the new code in situations involving both small and large uniform conductivities. 

The ID tests of the propagation of a large amplitude circularly-polarized Alfven wave show 
the new resistive relativistic MHD code reproduces ideal relativistic MHD solutions when the 
conductivity a is high and that the code has 2 nd order accuracy. The ID self-similar current sheet 
tests indicate that the analytical solution in the moderate conductivity regime is well described 
by the new code. In a simple MHD version of the Brio and Wu ID shock-tube test, the code is 
stable using Strang's splitting technique even when the conductivity is high (<r ~ 10 6 ). In the 2D 
cylindrical explosion tests, at the high conductivity a > 10 3 , the results recover the solution from 
ideal RMHD. The results of the Kelvin-Helmholtz instability test show that the growth rate of KHI 
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Fig. 16. — Evolution of (a) the maximum of the x-component of the velocity (v x ), (b) the volume- 
averaged magnetic energy (B 2 ), and (c) the reconnection rate as a function of time for the 2D 
relativistic magnetic reconnection test using an ideal EoS with T = 4/3 (solid lines) and a TM EoS 
(dashed lines). 

is independent of the conductivity, except for very low conductivity (a — 0). However, magnetic 
field amplification via stretching of the main vortex developed by KHI strongly depends on the 
conductivity. The effect of conductivity on magnetic field amplification via KHI in 3D is a topic 
for future study. 

EoSs proposed by Mignone et al. (2005) and by Ryu et al. (2006), which closely approximate 
Synge's single-component perfect relativistic gas EoS, have been incorporated in the code. In 
the limit of non-relativistic and ultra-relativistic temperatures, the equivalent specific heat ratio 
associated with the EoSs that approximate Synge's EoS appropriately changes from the 5/3 to the 
4/3 limits. 

The numerical tests studied the effect of the EoS on shocks, blast waves, the Kelvin-Helmholtz 
instability, and relativistic magnetic reconnection. The results provide a useful guide for future more 
specific studies of each topic. The tests confirm the general result that large temperature gradients 
cannot be properly described by an ideal EoS with a constant specific heat ratio. The results using 
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a more realistic EoS, which we have studied here, show considerable dynamical differences. The 
ID shock tube tests (Balsara Test2) show that the results obtained from the TM EoS and RC 
EoS cases are considerably different from the constant 7-law EoS case with 7 = 5/3. In the 2D 
Kelvin-Helmholtz instability tests, the non-linear behavior depended on the EoS, even though the 
growth rate of the KHI was almost the same. In reconnection tests, the approximate EoS cases 
resulted in a faster reconnection outflow speed and a larger reconnection rate than the ideal EoS 
case. We conclude that any studies of shocks, instabilities, and relativistic magnetic reconnection 
should use a realistic approximation to Synge's EoS. 
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